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Peaks in two-dimensional weak lensing (WL) maps contain significant cosmological information, 
complementary to the WL power spectrum. This has recently been demonstrated using N-body sim- 
ulations which neglect baryonic effects. Here we employ ray-tracing N-body simulations in which 
we manually steepen the density profile of each dark matter halo, mimicking the cooling and con- 
centration of baryons into dark matter potential wells. We find, in agreement with previous works, 
that this causes a significant increase in the amplitude of the WL power spectrum on small scales 
(spherical harmonic index t J> 1,000). We then study the impact of the halo concentration increase 
on the peak counts, and find the following, (i) Low peaks (with convergence 0.02 <CK pe ak <J0.08), re- 
main nearly unaffected. These peaks are created by a constellation of several halos with low masses 
(~ 10 12 — 1O 13 M0) and large angular offsets from the peak center ( J>0.5i? v ir); as a result, they 
are insensitive to the central halo density profiles. These peaks contain most of the cosmological 
information, and thus provide an unusually sensitive and unbiased probe, (ii) The number of high 
peaks (with convergence fv pe ak ;>0.08) is increased. However, when the baryon effects are neglected 
in cosmological parameter estimation, then the high peaks lead to a modest bias, comparable to that 
from the power spectrum on relatively large-scales [i < 2000), and much smaller than the bias from 
the power spectrum on smaller scales (£ > 2,000). (hi) In the 3D parameter space (a"g,f2 m ,w), the 
biases from the high peaks and the power spectra are in different directions. This suggests the possi- 
bility of "self-calibration" : the combination of peak counts and power spectrum can simultaneously 
constrain baryonic physics and cosmological parameters. 

PACS numbers: PACS codes: 98.80.-k, 95.36. +x, 98.65. Cw, 95.80,+p 



I. INTRODUCTION 

Weak gravitational lensing (WL) by large-scale cos- 
mic structures has emerged as one of the most promising 
methods to constrain the parameters of both dark energy 
(DE) and dark matter (DM) (e.g. ref. [T]; see also recent 
reviews in refs. [21 13])- WL was first detected over two 
decades ago [3]. It has recently matured to deliver cos- 
mological constraints; in particular, the COSMOS survey 
has provided independent evidence of the accelerated ex- 
pansion of the Universe [5J [BJ. Over the next decade, 
revolutionary large WL datasets are expected to be avail- 
able. The LSST survey will cover 20,000 square degrees 
with multi-band imaging suitable for weak lensing, and 
other large surveys will cover several thousand square 
degrees [7]. These WL datasets will be a rich source of 
cosmological information, going beyond traditional two- 
point statistics such as the power spectrum. 

Lensing peaks, defined as local maxima in two- 
dimensional WL maps, provide statistics that will be au- 
tomatically available and particularly straightforward to 
measure in lensing datasets. Since they probe the under- 
lying 3D density fluctuations on small, non-linear scales, 
they are sensitive to non-Gaussian features. As a cos- 
mological probe, the peak counts are therefore comple- 
mentary to the WL power spectrum, and are similar to 



galaxy cluster counts. A major advantage of peaks and 
a motivation for their use is that they avoid the issue 
of having to identify genuine bound clusters and mea- 
sure their masses. Accurate modeling of peak statistics 
requires large numerical simulations, which are now be- 
coming feasible. 

In the past few years, interest in lensing peaks and 
other, closely related statistics has increased signifi- 
cantly. 1 The probability distribution function of the con- 
vergence [9], and its cumulative version, the fractional 
area of "hot spots" on convergence maps [TU] are sim- 
ilar to peak counts in the high- convergence limit and 
have been shown to have useful cosmology sensitivity. 
The latter statistic is also known as Vo, one of the three 
Minkowski functionals (MFs) for two dimensional thresh- 
olded fields. MFs are related to peaks and had been pro- 
posed as a weak lensing statistic pTj [12] . More recently, 
ref. [13] constructed an analytical approximation to V2, 
the genus statistic (which also corresponds to peak counts 
in the high-threshold limit). The full set of Minkowski 



To our knowledge, lensing peaks were first considered as a cos- 
mology probe in the early ray-tracing simulations by ref. [8], 
which studied the f2 m — dependence of the peak counts. 
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functionals in the context of weak lensing was studied ex- 
tensively both theoretically [14 and in ray-tracing sim- 
ulations [H]. Finally, peak counts have also been stud- 
ied in wavelet space [TB], and found to break the usual 
(erg, f2 m )-degeneracy that exists between models from the 
power spectrum alone. 

Preliminary studies [17l [18] that defined peaks as lo- 
cal maxima were based on 2D projections of the 3D 
mass distribution in low-resolution simulations. WL peak 
counts with ray-tracing were subsequently studied by 
refs. [TH [20] and more recently in ref. [21]. These were 
based on simulations with better mass resolution, and 
revealed that low-amplitude peaks (which typically do 
not correspond to single collapsed dark matter halos) 
contain most of the cosmological information. Various 
other aspects of WL peak counts have been further ex- 
plored. Peaks have been shown [551 H3] to constrain the 
primordial non-Gaussianity parameter /nl- Ref. |24j in- 
vestigated the origin of the cosmologically important low 
peaks, and found that they are typically caused by a 
constellation of 4-8 low-mass halos. Ref. [TS] and [25] 
demonstrated that cosmological constraints from peaks 
can be tightened by combining several angular smooth- 
ing scales. WL peak counts have also been directly com- 
pared and found superior to two other commonly used 
non-Gaussian statistics, skewness and kurtosis |26j . Fi- 
nally, [27] study the effect of masked regions on shear 
peak counts and show that using Karhunen-Loeve analy- 
sis can mitigate biases on peak count distributions caused 
by masks, and that it can reduce the number of noise 
peaks. 

A common limitation to all of the above works is that 
they are based on cold dark matter (CDM) simulations 
(where simulations have been used), neglecting baryons 
and all astrophysical processes. This leaves an incom- 
plete description of the potential fluctuations and the 
corresponding lensing signatures. This is a particular 
concern since previous work has shown that the cooling 
and condensation of baryons inside dark matter halos 
can change the total matter distribution and has a large 
impact on the WL power spectrum on small scales (e.g. 
refs. [2"8Tl30| ). Furthermore, in astrophysical models that 
include feedback from active galactic nuclei (AGN) and 
supernovae, in addition to cooling and star formation, 
the matter distribution can be affected well outside dark 
matter halos, modifying the 3D matter power spectrum 
[3"T] . as well as the 2D WL power spectrum [32J out to 
large scales. 

Motivated by these findings, here we quantify the ef- 
fect of baryons on the statistics of WL peaks. Realistic 
modeling of the astrophysics, using hydrodynamical sim- 
ulations, remains challenging, both in terms of including 
all of the relevant physical processes correctly, and also in 
terms of computational scale. However, previous studies 
have shown that the cooling and condensation of baryons, 
and the resulting impact on the total (gas+DM) density 
profiles of halos, can be modeled by simple modifications 
to the halo density profile [2M501 1321 [53] . In particular, 



ref. |33| finds that a simple increase in the concentra- 
tion parameter of the universal NFW [31] profile can be 
a good approximation to the results of hydrodynamical 
simulations, and can account for the changes in the WL 
power spectrum |30) . We therefore follow this prescrip- 
tion, and manually steepen the density profile of each 
individual DM halo identified in our N-body simulations. 

This method is simple to use, and allows us to quantify 
the effects of gas cooling. A similar approach could be 
followed to model the effect of AGN feedback and other 
processes, but we leave this to future work. In this paper, 
we make a large change to the concentration (increasing 
it by 50%, compared to the 36% increase that was found 
to match simulation results 30J ) , thus intentionally am- 
plifying the impact of baryon cooling on the weak lensing 
statistics. 

The main goal of this paper is to investigate the bias 
in the cosmological parameters w, eg and Q m when peak 
counts and power spectra are fit neglecting baryonic ef- 
fects. Our results suggest that the bias from the peaks 
is lower, and also in a different direction, compared to 
the power spectrum. This suggests that the power spec- 
trum and peak counts can be combined to "self-calibrate" 
WL surveys, i.e. to fit cosmological parameters simulta- 
neously with parameters describing the halo profiles. 

The rest of this paper is organized as follows. In § [TTJ 
we describe our calculational procedures, including the 
creation of the WL maps through ray-tracing in N-body 
simulations, identifying halos and modifying their density 
profiles, and creating "baryonic" versions of WL maps. 
This section also presents our statistical methodology to 
compare maps, and our Monte Carlo procedure of es- 
timating confidence contours and the biases caused by 
neglecting the baryonic effects. In § |III[ we present our 
results, which include the effect of baryon cooling on the 
peak counts and on the power spectrum, as well as the 
biases caused by neglecting these effects when fitting for 
the three cosmological parameters w, ag and Q m . In 



IV we offer a detailed discussion of our main results, 



as well as of several caveats and possible extensions. Fi- 
nally, in § [V] we summarize our main conclusions and the 
implications of this work. 



II. METHODOLOGY 
A. N-body simulations and WL maps 

The cosmological N-body simulations of large-scale 
structures and ray-traced weak lensing maps used in this 
paper are the same as those in our earlier work [151 124j . 
We refer the reader to these publications for a full de- 
scription of our methodology; here we review the main 
features and describe the new features we have imple- 
mented to model the baryonic effects. We intend to make 
our data products publicly available in the future [35ll36j . 

A total of 80 CDM-only N-body runs were made 
with the Inspector Gadget lensing simulation pipeline 
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[551 155] at the New York Blue supercomputer. NY Blue 
is part of the New York Center for Computational Sci- 
ences at Brookhaven National Laboratory /Stony Brook 
University. 

Our suite of 7 cosmological models includes a fiducial 
model with parameters {£l m = 0.26, fl\ = 0.74, w — 
-1.0, n s = 0.96, cr 8 = 0.798, H = 0.72}, as well as 
six other models. In each of these six models, we varied 
one parameter at a time, keeping all other parameters 
fixed at their fiducial values; we thus have WL maps in 
variants of our fiducial cosmology with w — {—0.8, —1.2}, 
cr 8 = {0.75,0.85}, and O m = {0.23,0.29}. Note that 
in the last case, we set J7a = {0.77,0.71} to keep the 
universe spatially flat. 

To produce the N-body simulations, we first created 
linear matter power spectra for the seven different cos- 
mological models with CAMB [37] for z — 0, and scaled 
them back to the starting redshift of our N-body sim- 
ulations at z = 100 following the linear growth factor. 
Using these power spectra to create initial particle posi- 
tions, the N-body simulations were run with a modified 
version of the public N-body code Gadget-2 [35] and its 
accompanying initial conditions generator N-GenIC. We 
modified both codes to allow the dark energy equation of 
state parameter to differ from its ACDM value (w ^ 1), 
as well as to compute WL-related quantities, such as co- 
moving distances to the observer, at each simulation cube 
output. Each simulation contains 512 3 CDM particles in 
a cubic box with a side length of 240/i _1 comoving Mpc, 
allowing a mass resolution of 7.4 x 10 9 ft. _1 M Q . Each of 
these runs took approximately 1.75 real clock days, using 
64 processors on the Blue Gene. 

In each of the six non-fiducial cosmological models, we 
ran 5 strictly independent N-body simulations (i.e. each 
with a different realization of the initial conditions) . To 
minimize the differences between two cosmologies arising 
from different random realizations, the initial conditions 
for each of those five simulations were matched across the 
cosmologies quasi-identically. This entails recycling the 
same random number when drawing mass density modes 
from the power spectrum for each cosmology (note that 
the power spectra themselves of course differ across the 
cosmologies) . In the fiducial cosmology, we ran 50 strictly 
independent simulations - the first set of 5 to match the 
other cosmologies quasi-identically as mentioned above, 
and an additional set of 45 to improve the statistical ac- 
curacy of the predictions in the fiducial cosmology (espe- 
cially the covariance matrices). For the "baryonic" maps, 
as described below, only the first of these 50 simulations 
was used. Table[l]lists the sets of simulations in our suite, 
along with their cosmological parameters, the number of 
independent N-body runs, and the number of pseudo- 
independent 12 deg 2 maps (see below). 

To create the raw WL convergence maps, we used a 
standard two-dimensional, flat-sky ray-tracing algorithm 
closely following [39] with minor modifications as de- 
scribed in detail in [211] ■ Earlier work with similar al- 
gorithms includes [4"DT(4"2"] . Cubes with particle posi- 
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0.23 
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1000 
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0.798 


-1.2 


0.26 


0.74 
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1000 


w = -0.8 


0.798 


-0.8 


0.26 


0.74 
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1000 


erg = 0.75 


0.750 


-1.0 


0.26 


0.74 


5 


1000 


as = 0.85 


0.850 


-1.0 


0.26 


0.74 


5 


1000 


NFW-replaced 


0.798 


-1.0 


0.26 


0.74 


1 


200 


"baryonic" 


0.798 


-1.0 


0.26 


0.74 


1 


200 



TABLE I: Our weak lensing map sets, including the cosmolog- 
ical parameters, the number of underlying independent N-body 
simulations, and the number of weak lensing maps in each set. 



tions from the N-body simulations were output every 
80/i _1 Mpc in the radial (redshift) direction. While sev- 
eral independent simulations were used to make each 
map, boxes from the same simulations had to be recy- 
cled multiple times. The data cubes at each redshift 
were therefore randomly shifted, sliced, and rotated (by 
multiples of 90 degrees) to produce pseudo-independent 
realizations. The particles were projected onto two- 
dimensional density planes perpendicular to the central 
line of sight of the map. The triangular shaped cloud 
(TSC) scheme [33] was used to place the particles on a 
grid on these 2D density planes. The Poisson equation 
was then solved in Fourier space to convert the surface 
density into a gravitational potential. The deflection an- 
gles and convergences were calculated at each plane for 
each light ray from the first and second transverse spatial 
derivative of the 2D gravitational potential, respectively. 
Between density planes, the light rays were assumed to 
travel in straight lines; 2048 x 2048 light rays were fol- 
lowed in this fashion for each convergence map. A total 
of 1,000 pseudo-independent, 12 deg 2 convergence maps 
were produced in each CDM-only cosmological model; 
200 maps were produced in the sets used to study the 
systematic effect of baryons (see below). The number 
of maps in each map set is given in the last column of 
Table [H 

To create the final simulated WL maps, for simplicity, 
we assumed that the source galaxies are confined to the 
single redshift z s = 2. Ellipticity noise from the random 
orientations of the source galaxies was added to the con- 
vergence maps pixel by pixel, drawing from a Gaussian 
distribution corresponding a conservative source galaxy 
surface density of n ga i — 15 galaxies/ arcmin 2 , and an 
r.m.s. noise in one component of the shear of <7 7 = 0.22. 
Once noise was added, we smoothed the maps with a fi- 
nite version of a 9q — larcmin 2D Gaussian filter. This 
corresponds to the single most informative and smallest 
angular scale on which we trust our maps, based on com- 
parisons with the power spectrum from theoretical pre- 
dictions (see [15] for details). Combining several smooth- 
ing scales would tighten the overall constraints and could 
change the biases. We expect that this would not change 
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the qualitative results in this paper, but we leave an in- 
vestigation of this issue to future work. 

Peaks arc defined as local maxima on the pixelized final 
mock WL maps, and are counted in a straightforward 
fashion, with the convergence in the central pixel K poa k 
identified as the "height" of each peak. Power spectra 
are measured numerically from the Fourier transforms of 
the same maps following standard techniques. 

The differences between cosmological models, and 
thus the cosmological parameter dependence of the peak 
counts and the power spectra, are computed using pairs 
of the CDM-only 5-simulation map sets. Additionally, 
the covariance matrices of both observables - peak counts 
binned by their height, and power spectrum in bins of I - 
were computed from the 45-simulation fiducial map set; 
we have found that this was necessary for better accuracy 
(see [TS] for details). Note that we do not consider the 
dependence of the covariance matrices themselves on cos- 
mology (in principle, this dependence could help improve 
constraints). 

Finally, in order to study the impact of baryons, we 
create two more sets of WL maps. These are based on 
the fiducial CDM model, and a single N-body simula- 
tion. In this simulation, we identify all the DM halos in 
the 3D simulation cubes, and replace them with spher- 
ically symmetric halos with analytic NFW [33] density 
profiles. The ray-tracing procedure is then repeated ex- 
actly as before, and a new set of 200 "NFW-replaced" 
maps is created. This set is used to predict the expecta- 
tion values of the peak counts and the power spectra in 
the fiducial cosmology, in the absence of baryons. Finally, 
beginning with the same 3D data, the concentration pa- 
rameter "c" of each spherical NFW halo is increased by 
50%, and the ray-tracing and map-making procedure is 
once again repeated, to create a corresponding set of 200 
"baryonic" maps. This last set is used to find the expec- 
tation values of the peak counts and the power spectra in 
the fiducial cosmology, in the presence of baryon cooling. 
This procedure guarantees that any differences between 
an individual "NFW-replaced" map and the correspond- 
ing "baryonic" map is caused only by the changes in the 
halo concentration (and not from the halo replacements 
or from having different random realizations). Note that 
ideally we would want 1,000 pseudo-independent realiza- 
tions of the baryonic maps, to match the number we have 
for the CDM-only map sets. However, to keep computa- 
tional costs feasible, we have produced only 200 baryonic 
maps from a single N-body run. Nevertheless, we repli- 
cate each of these 200 baryonic maps 5 times, each time 
adding a different random noise realization. This pro- 
vides a better statistical sampling of the noise (in partic- 
ular, a more accurate determination of the average effect 
of the noise). 



B. Identifying and replacing dark matter halos 

In this section we describe in detail how we identify, 
remove, and re-insert halos into the N-body simulations. 

Once the N-body simulations are generated, we use 
the publicly available AMIGA halo finder ( 44 , hereafter 
AHF) to identify collapsed halos in our N-body runs. 
Its output consists of the 3D positions and the tagged 
set of particles belonging to each halo, the central posi- 
tion of the halo (defined as the local maxima of the den- 
sity field), as well as the total number of halo particles 
inside spherical shells at several different discrete radii, 
from which the spherically-averaged density profiles can 
be calculated. Depending on the redshift, AHF identifies 
a total of 2.5 x 10 5 — 2.6 x 10 5 halos in our box, with 
masses in the range 1.5 x 10 11 h~ 1 M e - 5 x 10 14 h~ 1 M Q . 
We then remove all the particles belonging to all of the 
halos from the N-body simulations, and add back the 
density profiles fitted to the identified halos, assuming 
halos are spherically symmetric and described by NFW 
profiles. The details of the fitting will be given in §11 C| 
below. 

In the procedure of replacing halos, three non-trivial 
points regarding mass conservation need to be clarified. 
The first point concerns sub-halos. When a parent halo 
and a sub-halo share a common structure, the halo finder 
saves the shared particles in both the parent and the sub- 
halo profiles. As a result, in our procedure, the particles 
within halos are removed one time, but the shared parts 
of halos are added back repeatedly. We have found that 
this can cause an artificial 5% increase in total (halo plus 
subhalo) mass. To avoid this problem, and to conserve 
mass, we sort the halo catalog by mass. Beginning with 
the lowest-mass halo, we consider halos in this ranked list 
one-by-one, always comparing the position of the center 
of a halo with the positions of all of the halos further 
down the list (with higher masses) . If the center of a halo 
is found to fall inside the virial radius of any of the higher- 
mass halos, we subtract the mass of the smaller halo from 
that of the larger one. When we re-insert the larger halo 
into the 3D simulation box, we use the analytical NFW 
profile with the reduced mass. 

The second point concerns discretization of the ana- 
lytic NFW profiles. In our map-making procedure, the 
projected halo density profiles are evaluated only at the 
discrete grid points on our 2D density plane. As a re- 
sult, each halo effectively contributes a total mass given 
by a "2D trapezoidal integral", rather than its actual 
(analytically calculable) mass. In principle, the actual 
halo profiles are known, and the projected mass could 
be resolved to arbitrary accuracy; in practice, we have 
found this to be computationally too expensive. Instead, 
in order to conserve the total mass, we normalize the 
discretized surface density profile by multiplying by the 
ratio of the actual mass of the halo and the total pro- 
jected mass of the discretized profile. 

Finally, a third point concerns halos that are "trun- 
cated", either in the transverse direction (because they 
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are too close to the edge of the 12 deg 2 map) or in the 
rcdshift direction (because they arc too close to the back 
or front side of the underlying 3D simulation slice; note 
that the raw 3D simulation cubes have periodic bound- 
ary conditions, but the slices used to make the maps do 
not [H]). In the former case, we add the projected halo 
density profiles as we do for the normal halos, apply- 
ing periodic boundary conditions. In the latter case, the 
surface density of a truncated halo is taken by summing 
the discretized density profile over the region inside the 
simulation box. We do not re-normalize the profiles of 
these truncated halos in the redshift direction, because 
it is difficult to determine how much mass was actually 
lost in the truncation. However, we have shown in our 
previous study [24] that the effect of edge halos can be 
safely ignored for studying the peak counts. 

After accounting for the sub-halos and normalizing the 
discretized profiles, the fractional difference between the 
total mass removed and the total mass added back is less 
than 0.2%. 

To check the ultimate accuracy of the halo replacement 
procedure, in Figure [T] we plot the fractional difference 
in the peak counts and the power spectrum caused by 
the halo replacement in the fiducial cosmology. The data 
points are averaged over 100 realizations. For this check, 
we used the noiseless versions of the WL maps. As the 
upper panel shows, the halo replacement works very ac- 
curately for peaks below k % 0.06. The difference is still 
smaller than 10% for peaks below n < 0.1, however there 
is a systematic decrease in the number of peaks that be- 
comes significant for higher amplitude peaks. This de- 
crease is caused by our treatment of the sub-halos. High 
peaks are typically caused by a single massive halo. Mas- 
sive halos contain many sub-halos, and because of the 
substantial mass loss due to tidal stripping, sub-halos 
are found preferentially away from the center of the host 
halos [45 . Recall that when we add back an NFW par- 
ent halo, in order to conserve mass, we subtract the total 
mass in subhalos. When we perform this subtraction, 
we assume, for simplicity, that the mass in the subhalos 
is distributed radially in a smooth fashion, following the 
density profile of the parent halo. This procedure there- 
fore over-subtracts mass near the halo center, and under- 
subtracts mass near the halo outskirts. The overall effect 
of our procedure is to decrease the central density (by a 
few percent). Since the line of sight to a high peak typi- 
cally passes very close to the center of the corresponding 
massive halo, this diminishes the height of the peak, and 
results in fewer high peaks. A similar explanation holds 
in the case of the power spectrum, shown in the lower 
panel. The halo replacement works very accurately for 
large scales, with AP/P < 2% for I < 1000. But at 
I >5,000, where the one-halo term dominates the power 
spectrum, there is a significant (~ 6%) decrease in power. 

We emphasize that the inaccuracies due to halo re- 
placement are quantified here only as a reassurance that 
we did not, in the process, gravely modify the WL observ- 
ables. We are only interested in the impact of baryons; 



indeed, the impact of baryon cooling in the real 3D halos 
should be similar to the increase in concentration param- 
eters for the spherically symmetric halos. 
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FIG. 1: The effect of replacing halos in our N-body simu- 
lations by spherically symmetric NFW halos. In the upper 
panel, we show the fractional difference (AN) / (N) in the 
peak counts caused by this halo replacement, in convergence 
bins of width An — 0.00675. In the lower panel, we show 
the corresponding fractional difference in the power spectrum 
(AP K ) / (P K ) , in logarithmic bins of the spherical harmonic in- 
dex with a width A\og£ = 0.2. Data points are averaged over 
100 realizations. The error bars in the top panel are estimated 
as the standard deviation of AN , divided by (N) ( and simi- 
larly as the standard deviation of AP K , divided by (P K ), for 
the power spectrum). The source galaxies are assumed to be at 
redshift z s = 2, and the convergence maps have been smoothed 
with a 1 arcmin Gaussian filter. 



C. Fitting halo density profiles 

In this section, we describe our procedure for the fit- 
ting the halo density profiles and quantify how well this 
procedure works. 

We assume all halos and subhalos follow a universal 
NFW profile [34], 

Ps 



/Onfw(r) 



(r/r s )[l + (r/r s )Y< 



(1) 
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where r is the radius from the halo center, and r s and 
p 8 are a characteristic radius and density. The profile is 
truncated at R\soi inside which the mean overdensity is 
180 times the matter density of the universe at redshift 
z. This convention is consistent with the settings in the 
AHF halo finder. The concentration parameter is given 
by ciso = Riso/r s . 

The NFW density profile is uniquely defined by the 
concentration parameter ciso, and the normalization fac- 
tor p s . To obtain Ciso, we chose to fit the normalized 
cumulative density profile F(x), i.e. the fraction of mass 
within a certain radius r, as a function of x = r/Riso, 

p, x j _ log(ci 80 x + 1) - Cwqx[ (1 + CisqO:) ^ 
log(ci80 + 1) - Ciso/(l + Ciso) 

Here Ri$o is taken as the actual radius of halos found by 
the halo finder. We fit F(x) to find Ciso and calculate 
the normalization factor p s as: 



M 



Ps 



LSI) 



47rrf (log(ci 80 + 1) - ciso/ (1 + ciso)) ' 



(3) 



Here Miso is the actual halo mass returned by the halo 
finder. For some halos, we find that c 18 o does not con- 
verge to an acceptable number (1.1 < ciso < 50), or the 
halo mass is too low and the halo finder cannot compute 
a radial profile. We distinguish these as unfitted halos, 
and for these halos we use an analytical formula to obtain 
Ciso, adapted from Eq. (5) in ref. |46| : 

c 200 = 4.67 (Af 200 /10 14 / l - 1 M Q )- all /(l + z). (4) 

Eq. Q adopts a different convention from ours - the 
profile is truncated at the radius inside which the mean 
overdensity is 200 times the critical density of the uni- 
verse. We therefore extend the NFW profile for a halo 
with mass M200 outward to a radius where the mean inte- 
rior density is 180 times the background matter density. 
This extrapolation yields a relation between the input 
A/iso and M200 that depends on ciso; using Eq. Q and 
the definitions of the various quantities above, the pa- 
rameter ciso = Ri8o/ r s can be found iteratively. 

The unfitted halos account for ~ 2% of all halos in 
our simulations. This is unsurprising, since, for example, 
halos that are undergoing major mergers, or have not yet 
relaxed from a recent major merger, have no reason to 
follow NFW profiles [46]. 

We find that the rest of our halos follow the NFW 
profiles accurately. In Figure [2j we show the spherically 
averaged density profile p(r)r^, for a halo with a mass 
of 8.7 x 10 13 h~ 1 M & at redshift z = 0, along with the 
corresponding best-fit NFW profile. The quality of the fit 
shown in this figure is typical of halos in our simulations. 

As a quantitative test of the accuracy of our fitting, we 
define a statistic 



T = 



1 n 



F, - F( Xi 



F 



(5) 



Here n is the number of data points (i.e. in radius) for 
one halo, Fi and F(xi) is the «th data point and the 
fitted value at that radius. We plot the distribution of T 
in Figure [3j The figure demonstrates that typically, the 
NFW profiles are accurate to within 5%, although there is 
a tail of outliers extending to poorer fits. As mentioned 
above, approximately 2% of our simulated halos could 
not be fit by NFW profiles at all (and are not shown in 
this figure). 
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FIG. 2: The figure shows the spherically averaged density pro- 
file p(r)r 3 (red crosses), for a halo with a mass of 8.7 x 
1O 13 /i -1 M0 at redshift z = 0, along with the corresponding 
best-fit NFW profile (dashed blue curve). The quality of the 
fit shown in this figure is typical of halos in our simulations. 
Here p m is the mean matter density of the universe and Riso 
is (approximately) the virial radius of the halo. The best-fit 
NFW profile was obtained by a least-squares fit to the normal- 
ized cumulative density profile (Eq. Ufy). 



D. Statistics 

We refer to the different statistics one can obtain from 
a 2D WL map — e.g. power spectrum, peak counts, etc. — 
as A/j. The index i in the components 2Vj of the vec- 
tor N labels different heights for the peaks, and differ- 
ent multipoles for the power spectrum. We divide the 
peak counts into 30 bins in height K pea k in the range 
0.023 < Kpcak < 0.08 and, we use a second set of 
30 bins above K pea k > 0.08. Similarly, we divide the 
power spectrum into 30 multipole bins t in the range 
100 < t < 2 x 10 4 . For computational reasons, the power 
spectrum is first pre-computed and stored for 1,000 "pre- 
bins" spaced linearly in multipole 100 < £ < 1 x 10 5 and 
is only later combined into the 30 larger bins used for the 
final computation. 



This discussion closely follows ref. I15| . 
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FIG. 3: The distribution of the statistic T, defined as the frac- 
tional deviation between the fitted and the actual normalized 
cumulative density profile (Eq. ^3^), averaged over radii, for 
the halos in our N-body simulations. The figure demonstrates 
that typically, the NFW profiles are accurate to within 5%, 
although there is a tail of outliers extending to poorer fits. 
Approximately 2% of our simulated halos could not be fit by 
NFW profiles at all (and are not shown in this figure). 

To constrain cosmology, we are interested in the true 
ensemble average, over an infinite number of realiza- 
tions within a single cosmology (hereafter denoted by 
brackets ( )) as a function of cosmological parameters 
p = {fi m , w, erg}, as well as the ensemble covariance. 3 
Since these are not available, we estimate them from 
a finite number of our simulations. Averaging over the 
pseudo-independent map realizations within a given cos- 
mology, we can obtain an estimate for the ensemble av- 
erage by 

_ i H 

(jv,(p))wjv i ( P )= s x; jv *( r 'P)' ( 6 ) 

r— 1 

where A^(r, p) is the statistics vector measured in a sin- 
gle map and r runs over our R = 1,000 maps. We call 
this estimate the simulation mean. It differs from the 
true ensemble average both because of the limited num- 
ber of realizations and also because of the limitations 
inherent in our simulations, such as limited resolution. 
In the absence of a fitting formula for the peak counts in 
the non-Gaussian case (analogous to the power spectrum 
formula from [47] in the nonlinear regime) the simulation 
mean serves as our proxy for theoretically predicted peak 
counts. 4 



3 More generally, one could utilize the full probability distribution 
of a statistic, not just its average, as well as the cosmology- 
dependence of all higher-order correlations, not just the covari- 
ance; we will not investigate these issues in the present paper. 

4 The ensemble average for the peak counts can be computed ex- 
actly for a Gaussian random field |48l ; however, this is a poor 
approximation to the lensing peaks 1241 . 



We can form this estimate only for the 7 selected cos- 
mologies where we have run simulations (Table [I]). For 
cosmologies with other parameter combinations p, we 
have to interpolate (and in a few cases extrapolate) be- 
tween these 7 points in the 3D parameter space. We 
thus construct a first-order Taylor expansion around our 
fiducial cosmology: 

F,(p) « F«(p< ' NFW )) + 

+ (a)_ (0) (Pc*-p ( a } )-(7) 

a Pa Pa 

The index a — 1,2,3 refers to either Cl m ,w, or erg, and 
the sum counts through all 3 parameters. The fraction 
on the right-hand-side of Eq. Q is a finite difference 
derivative along the direction of the parameter a, com- 
puted using the fiducial cosmology (p^) and a cosmol- 
ogy in which the parameter p a was changed from the fidu- 
cial case (p^). As mentioned above, the 5-simulation 
CDM map sets with quasi-identical initial conditions 
were paired and used for the computation of this deriva- 
tive, but the simulation mean for the NFW-replaced set 
was used as the expansion point of the Taylor series (the 
first term on the RHS). We indicate this explicitly in 
Eq. Q by the superscripts "CDM" and "NFW" . 

If this non-fiducial cosmology is chosen such that 
Pa^ — p^ is positive for all three parameters, we call 
it a "forward derivative" , if it is negative, we call it a 
"backward derivative" . We will also refer below to a "2- 
sided derivative" , which switches automatically between 
the forward and backward cases as needed for each pa- 
rameter, corresponding to the octant in the 3D parameter 
space where the interpolation is performed. For simplic- 
ity, we use the backward derivative throughout this paper 
unless explicitly noted otherwise (and below we highlight 
some important differences in our results obtained with 
other derivatives). 

Similarly to the simulation mean, we estimate 
the ensemble covariance matrix from the simulations, 
Cov(Ni,Nj) fa Cij, where 



Cij (P) = £ l N * ^ P) - N * (P)] M ^ P) - N i (P)l ■ 

r— 1 

This covariance matrix contains contributions both from 
the sample variance of the signal and from the galaxy 
shape noise. In the case of the power spectrum, the 
Gaussian random noise in different I bins should be un- 
corrected, and therefore contribute only diagonal terms 
(in practice, the cross-terms are small but nonzero in our 
finite set of realizations). However, adding noise to the 
maps has a nonlinear effect on peak counts, and intro- 
duces off-diagonal terms, as well. As mentioned above, in 
practice, we only evaluate and use the covariance matrix 
in our fiducial 45-simulation map set. 



E. Monte Carlo parameter estimation contours 

Each of our WL maps spans a 12 deg 2 field of view, 
yet we wish to obtain parameter contours for a 20,000 
deg 2 full-sky survey, such as LSST. Thus we employ a 
parametric bootstrapping approach to generate approxi- 
mations to full-sky maps. In this procedure, we randomly 
select a map from our set of 1,000 maps, with replace- 
ment, 20, 000/12 rj 1, 667 times. The effective solid angle 
of this larger "composite map" is then 20,000 deg 2 , as 
desired. The map, of course, is not a true composite - 
we merely compute the observables for each of the 1,667 
patches individually, and then add or average over them 
to get their values for the full-sky map. We create 10,000 
such full-sky maps to obtain smooth parameter contours 
in our Monte Carlo procedure. While this process can- 
not generate information not present in the 1,000- map 
set, the bootstrap procedure has a large benefit. Individ- 
ual 12 deg 2 maps constrain parameters only poorly, and 
large fluctuations in JVj require parameter extrapolations 
far outside the range of our simulations. The averaging 
in the bootstrap procedure suppresses these fluctuations, 
so the Ni are similar to those from the simulations. 

To estimate the cosmological parameter error contours 
from the set of baryonic maps, we use ^-minimization 
to find the best-fit cosmological parameter combination 
for each of the 10,000 bootstrapped full-sky maps. Here 
X 2 is defined by 



a modest extrapolation in the parameter space. In con- 
trast, the scatter of best- fit points from the original small 
12 deg 2 maps was so large that a significant fraction of 
the best-fits lie outside the simulated region, requiring 
large extrapolations. Bootstrapping also eliminates the 
need to re-scale the contours to correspond to a full-sky 
survey. 



III. RESULTS 
A. Effect of baryons on the statistics 




100 1000 10000 

1 



x 2 (r b ,p) = AJVi(r 6l p) [Cov" 1 ^)^ ANj(r b , p) 

(9) 

where 

ANi(r b ,p) = JV 4 (r 6j p<°: bBr » on )) - (^(p)>- (10) 

Note that r b here counts the full-sky maps and therefore 
runs over r b — 1, . . . , 10, 000. For each Monte Carlo real- 
ization r b , we minimize \ 2 with respect to p using a sim- 
ulated annealing algorithm. This is a popular technique 
based on Markov Chain Monte Carlo with a decreasing 
"temperature", which helps the algorithm to get out of 
local minima (we had success with this algorithm in our 
earlier work [T5]). We emphasize that the matrix Cov 
in Eq. Q can be reasonably arbitrary and does not even 
have to be the covariance matrix; this makes this method 
robust with respect to noise in the covariance matrix es- 
timate [T5] . 

The bootstrapping procedure explained above returns 
10,000 sets of best-fit parameters for each full-sky real- 
ization; the distribution of these best-fit points can be 
used to draw confidence contours at desired confidence 
levels (68% in this paper, unless stated otherwise) based 
on the density of the best-fit points. 

We indeed find that the best-fit points lie within the 
parameter range of our simulations, so that interpola- 
tion can be used during the best-fit procedure. The only 
exception is the high-£ power spectrum, which required 



FIG. 4: The fractional difference (AP K ) / (P K ) in the power 
spectrum caused by the 50% boost in the concentration pa- 
rameter of each NFW halo, shown in models without (red 
solid curve) and with galaxy shape noise (blue dashed curve; 
the noise corresponds to a source galaxy surface density of 
rigai = 15 arcmin -2 at redshift z s =2). The data points are 
averaged over 200 realizations in logarithmic bins of width 
A\og£ = 0.2. Error bars are estimated as the standard de- 
viation of the AP K , divided by (P K ) . We have artificially in- 
creased the error bars by a factor of 10 for visual clarity. The 
smallness of the error bars indicate that the impact of baryons 
on the power spectrum is highly systematic, i.e. very similar 
in each realization. 

We begin the presentation of our results with the effect 
of baryons on the statistics themselves. In Figure |4j we 
show the fractional difference (AP K )/(P K ) in the power 
spectrum, caused by the 50% boost in the concentration 
parameter of each NFW halo, shown in model with and 
without galaxy shape noise. The data points are aver- 
aged over 200 realizations in logarithmic bins of width 
Alog<? = 0.2. As the figure shows, in the noiseless case, 
a 50% increase in concentration causes a strong increase 
in the small scale power, by about 20% at £ = 10 4 . The 
effect is much smaller on large scales (I <2000). This 
agrees with the qualitative conclusions in ref. [3D], who 
used a similar toy model. The drop beyond £ >3 x 10 4 
is due to our 1-arcmin smoothing. In the noisy case, the 
fractional difference has a turnover at the much larger 
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FIG. 5: Similar to Figure [^J except we show the fractional 
difference in the peak counts, (AA r pea k}/(A ? peak}, caused by 
the same 50% boost in the halo concentration parameters. The 
data points are averaged over 200 realizations, and shown in 
convergence bins of width Ak = 0.00675. The error bars are 
shown at their original size. Compared to the power spectrum, 
the impact of baryons on the peak counts is less systematic, 
i.e. varies more significantly between realizations. 

scale of I « 3,000; this is the scale at which the galaxy 
shape noise starts to dominate. Note that the absolute 
difference (AP K ) in the power spectrum in the noisy case 
is exactly the same as in the noiseless case. Therefore 
even the small < 2% fractional difference in the noisy case 
can cause a strong bias in inferred cosmology when the 
baryonic power spectrum is fitted by the non-baryonic 
model (see detailed discussion below). 

In Figure [5j we show the fractional difference in the 
peak counts, ( AA poak )/ (-/Vp eak ) , caused by the same ar- 
tificial 50% boost in the halo concentration parameter. 
The data points are again averaged over our 200 maps, 
and shown in convergence bins of width Ak = 0.00675. 
As the figure shows, in both the noiseless and the noisy 
case, there is a strong increase in the number of high 
peaks, but very little change in the number of low peaks. 
The reasons why the low peaks are robust to the change 
in the concentration parameter will be discussed below. 



B. Cosmological constraints and biases 

We next present the biases in the inferred best-fit cos- 
mological parameters when the baryonic effects are ig- 
nored. We define the "low" and "high" peaks to have 
heights cr„oise < K pcak < 0.08 and K pcak > 0.08, where 
"noise = 0.023 is the r.m.s. of k from galaxy shape noise. 
We similarly define the power spectrum on "very large" , 
"large" , and "small" angular scales (or "very low" , "low" , 
and "high" I ranges) to be those with 100 < I < 1000, 
100 < I < 2000 and 2000 < t < 2 x 10 4 , respectively. 
Note that two of these ranges overlap partially. 

As described in § |II E[ we generate the distribution of 



best-fit points in the 3-dimensional parameter space of 
CTg, w, and fi TO . Projecting this 3D distribution in one of 
the 3 dimensions can then be used to define the 68% joint 
confidence contours on the remaining two parameters, 
marginalized over the third. As also described above, the 
baryonic model was compared to the non-baryonic NFW- 
replaced model, linearly interpolated with the backward 
derivative, using noisy maps with z s = 2 and 1 arcmin 
smoothing. We apply 30 nearly logarithmic bins to the 
whole range 100 < £ < 2 x 10 4 of the power spectrum. 
Note that this leaves fewer than 30 bins in each of the 
very low, low and high-£ ranges as defined above (namely, 
9, 17, and 13 bins, respectively) The nearly logarithmic 
bins are generated by choosing the bin boundaries so that 
the 1,000 pre-bins are assigned to 30 bins with nearly 
equal logarithmic spacing. For very low £, this results 
in bins that are linearly spaced, because there are not 
enough pre-bins to achieve truly logarithmic spacing. For 
the peaks, the bin boundaries are chosen so that each of 
the 30 bins for the low peaks contains the same num- 
ber of peaks, and likewise we use 30 equal-count bins for 
the high peaks. This split in the equal-count procedure 
avoids having too few bins for the high peaks. 

In Figure [6j we show the 68% confidence contours ob- 
tained from the peak counts, along with those from the 
power spectra, for a 2 x 10 4 -deg 2 full-sky survey. We have 
approximately 2,000 low peaks and 250 high peaks in ev- 
ery 12 deg 2 field. The offsets between the center of the 
underlying 3D distribution (i.e. the best-fit model in 3D) 
and the correct fiducial values are listed in Table [TTJ and 
correspond to the biases on the inferred cosmological pa- 
rameters. The contours in the 2D plots of Figure [6] are 
always marginalized over the third parameter that is not 
displayed, i.e. they are projections of the full 3D best-fit- 
point distribution onto the 2D plane, and the contours 
are calculated in this projected space. Recall that we 
only have 1000 (pseudo-)independent 12 deg 2 maps, from 
which 1,666 maps are drawn randomly with replacement 
in the bootstrapping procedure. As bootstrapping can- 
not improve the accuracy of the mean, it cannot improve 
the determination of the locations of the centers of the 
contours, and affects only their size. The "lcr" error of 
the centroids (in these 2D planes) should therefore be 
larger than the size of the contours (the 68% confidence 
levels from the LSST-like survey) by a factor of approx- 
imately ^1666/1000 « 1.3. This is in the limit that 
random noise dominates the dispersion in the best-fit lo- 
cations. In the opposite limit, i.e. if the dispersion is 
dominated by the variations among the underlying 200 
baryonic maps, then the peak likelihood locations in 2D 
would be larger by an additional factor of y/E. 

From Figure [6j we can roughly estimate that the joint 
2D errors from the peak counts extend over Sos « 0.005, 
Sw « 0.03 and SQ m « 0.005; the corresponding uncer- 
tainties in the locations of the best-fits in 2D are there- 
fore Sa' 8 w 0.0065, Sw' « 0.04 and SQ' n w 0.0065. A 
similar estimate can be done for the power spectrum. 
Given these uncertainties, we conclude that the low peaks 
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Observable 


Acs Aw Af2 m 


Low Peaks 


0.009 0.01 -0.001 


High Peaks 


0.061 0.10 -0.024 


Low-£ Pow. Spec. 


0.03 0.11 -0.01 


Very Low £ Pow. Spec. 


0.01 0.07 -0.00 


High-£ Pow. Spec. 


-0.23 -0.58 0.15 



TABLE II: The biases in the three cosmological parameters 
from the peak counts and power spectrum. The numbers in 
the table show the offset between the location of the best-fit 
in 3D and the correct fiducial cosmology along each of the 3 
parameters as, w, and fi m . 



and the very low-£ power spectrum are both robust to 
baryonic cooling, with the inferred biases consistent with 
zero; the exception is that low peaks have a small bias in 
erg. The high peaks and the high-£ power spectrum both 
have significant biases, Acr 8 ~ 0.061, Aw sa 0.10 and 
AO m k, —0.024 from peak counts, and Acrg ~ —0.23, 
Aw sa —0.58 and Af2 m w 0.15 from the power spectrum 
(see Table [II). The low bias, at least within the context 
of our mode , and the high information content of the 
low peaks suggests the value of including them in cosmo- 
logical analyses for dark energy (see further discussion in 



§ |IVB| below). 

Perhaps our most interesting result is that the biases 
from the peaks and the power spectra are in different di- 
rections in the 3D parameter space. The bias from the 
high peaks, in particular, is in a direction nearly oppo- 
site from the (much stronger) bias from the high-£ power 
spectrum; they are also in different directions from the 
\ow-£ power spectrum. This opens up the possibility 
for self-calibration, whereby baryonic parameters, such 
as the halo concentration parameter in our case, could 
be determined simultaneously with the cosmological pa- 
rameters. Only for the correct concentration parameter 
will the contours from the different observables align. 



C. Bias directions and goodness-of-fit 

In this section, we address two important questions 
related to the above results: (1) what is the quality of 
the biased best-fits? (2) what determines the direction 
of the biases? The first question is especially important, 
since in principle, a poor best-fit can reveal the presence 
of unaccounted-for parameters. 

In Figure [7j we show the peak counts and power spec- 
tra in the baryonic model, compared directly to the best- 
fit baryonless models. The four panels (top left, top right, 
bottom left, bottom right) show the results for low peaks, 
high peaks, low £ power spectrum and high i power spec- 
trum (following the definitions of the low/high ranges in 
§ |IIIB[ and using the corresponding best-fit points in Ta- 
blc |II[). In each case, the deviations from the baryonless 
fiducial cosmological model are shown. In each panel, the 
solid [red] curves show the effect of the baryons, and the 
dashed [blue] curves show the deviations that best mimic 



these curves, achieved in the best-fit baryonless models. 
As shown by the figure, in all four panels, we find an 
excellent fit, i.e. the best-fit curves agree well with the 
baryonic curves. Although the best fit somewhat over- 
predicts the number of peaks at the high-K end of the 
low peaks in the top left panel, we find a total x 2 = 30. 
Since we have 30 bins, the fit is good, with a reduced 
X 2 of unity (with somewhat better fits in the other three 
panels). This implies that when the baryonic effects are 
neglected when fitting the data, the cosmological param- 
eters can be biased, and this will not be obvious from the 
low quality of the fit. 

To investigate the origin of the bias directions, the 
long-dashed [green], dotted [pink] and dot-dashed [cyan] 
curves show the effect of changing a single one of the 
three parameters as, w, Q m at a time, to its best fit 
value, i.e., the bias in one parameter times the backward 
derivative of the observable with respect to that. (This 
means that the green, pink and cyan curves sum up to 
the blue curve.) 

In the case of the low peaks (top left panel), the best 
fit is driven almost entirely by the bias in cr$: the fig- 
ure shows that the change in as, by itself, can mimic 
the baryonic effects quite well. 5 Furthermore, the small 
residual would require decreasing the number of peaks 
at the high-K end (in bins 20 and higher), while keeping 
the peak counts in the lower-K bins unchanged. Such a 
change in the shape of the peak count distribution can 
not be accomplished by changing either w or f2 m ; as a 
result, these two other parameters remain essentially un- 
biased. 

The case of the high peaks (top right panel) is quite 
different, with degeneracies between all three parameters 
playing an important role in the bias. In particular, the 
changes in the counts due to as and fi m have a strong 
degeneracy, so that these two parameters work in the op- 
posite direction. 6 However, the degeneracy is not perfect, 
and leaves a nonzero residual - interestingly, this resid- 
ual can be mimicked essentially perfectly by a relatively 
modest change in w. 

Turning to the low-£ power spectrum (bottom left 
panel), the fit is again driven by the bias in cr 8 , but unlike 
in the case of low peaks, as alone cannot mimic the bary- 
onic effects. As a result, w and Q m play a role, working 
in the direction opposite to as, and producing a near- 
perfect fit. 

The most interesting case is the high-^ power spec- 
trum, shown in the bottom right panel in Figure [7| Here 
the degeneracy between as and fl m is so perfect that the 
changes in these parameters (with opposite sign) essen- 
tially cancel each other. As a result, the net change in 



5 Note that reducing the number of low peaks requires an increase 
in as; this somewhat counterintuitive result has been found [20| 
and explained in detail |24| in our earlier work. 

6 This erg — Sl m degeneracy in the abundance of ~ 3<r peaks has 
been noted earlier [19 24 . 
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FIG. 6: The joint 68% confidence contours on pairs of cosmological parameters (marginalized over the third), obtained from peak 
counts and the convergence power spectrum, in a full sky (2 x 10 4 deg 2 ) survey. The location of the correct fiducial cosmology 
is marked by a green dot in each panel. For reference, crosses mark the other cosmologies we simulated. 10,000 Monte-Carlo 
realizations of the baryonic model were fit by models neglecting the baryon effects, which can produce a bias in the inferred 
parameters. In the upper panels, results from the low (a n0 i se < K pca ^ < 0.08) and high (upc^ > 0.08) peaks are shown by 
solid red and solid blue contours. In the lower panels, results from the very-low (I < 1000,), low (I < 2000,), and high-l 
f2,ooo < i < 2 x io 4 ; power spectra are shown by dashed red, solid red and solid blue contours, respectively. The bias is 
strongest from the small-scale (high-i) power spectrum; it is nearly negligible for the low peaks and for the very-low £ power 
spectrum. 



the high-£ power spectrum can be attributed almost en- 
tirely to the change in w. This explains why there is a 
very large bias in w. The fact that the bias from the high 
peaks is driven by ogi whereas the bias from the high-i' 
power spectrum ends up being dominated by w, makes 
it less surprising that the biases from high peaks and the 
high i power spectrum are in very different directions. 



IV. DISCUSSION 

A. Robustness of the inferred bias 

How robust are the results to the number of bins? Our 
results about biases have not been optimized over the 
number and placement of bins. We have found, in gen- 
eral, that all the results listed in the last section hold 
qualitatively, as long as the number of bins used is not too 
small. In Table III we show the biases of the three cos- 
mological parameters as, w, Q m inferred from the peak 



counts, with three different numbers of bins, using noisy 
maps with z s = 2. Bin boundaries were chosen such that 
each bin contains the same number of peaks. The table 
shows that the biases from the low peaks are quite sta- 
ble; as the number of bins is varied from 20 to 30, the 



Observable 


bins 


A<T8 Aw AQ m 


Low Peaks 


20 
25 
30 


0.014 0.01 -0.004 
0.009 0.03 0.001 
0.009 0.01 -0.001 


High Peaks 


20 
25 
30 


0.057 0.16 -0.020 
0.069 0.07 -0.030 
0.061 0.10 -0.024 



TABLE III: The inferred biases of the three cosmological pa- 
rameters as, w, tt m from the peak counts, using three differ- 
ent number of bins. Low and high peaks are always defined 
to be those with heights between a n0 i se < K pca k < 0.08 and 
Kpcak > 0.08; with o n oise = 0.023. Bin boundaries were cho- 
sen such that each bin contains the same number of peaks. 



changes in the biases are always less than the uncertainty 
on the bias (5a' s « 0.0065, Sw' 0.04 and 5tt' m « 0.0065 
as mentioned above). High peaks have larger changes 
in the corresponding biases than the low peaks, but the 
sizes of these fluctuations are still comparable to the un- 
certainties. We conclude that our qualitative results from 
the peak counts are robust as long as >20 peak-height 
bins are used. 

How does the result depend on the finite difference 
derivatives? To compute the inferred biases, the bary- 
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FIG. 7: The figure compares the changes in the peak counts and the power spectra caused by the presence of baryons (solid 
red curves), and by the biased cosmological parameters in the best-fit baryonless models (dashed blue curves). For every curve, 
the corresponding values in the fiducial, baryonless model is subtracted, so that only the offsets from this fiducial model are 
shown. The four panels (top left, top right, bottom left, bottom right) show the results for low peaks, high peaks, low-l power- 
spectrum and high-i power spectrum, respectively. Bin numbers are shown on the bottom x axis, with corresponding k and 
t values indicated on the top of the figure (note that the binning is not linear in k and I; see § |77Ji3| of the text for a full 
description of the binning scheme). In general, the biased cosmologies are able to mimic the effect of baryons remarkably well. 
The long-dashed [green], dotted [pink] and dot-dashed [cyan] curves show the "decomposition" of the biased best-fit cosmology 
among the three cosmological parameters: they show the effect of changing a single one of the three parameters 
time to its best fit value. (The green, pink and cyan curves sum up to the blue curve.) 



onic models were fitted with the non-baryonic models, 
using linear interpolation and backward finite difference 
derivatives. Unless the model being evaluated is very 
close to one of the simulated cosmologies, there is no 
clear reason to prefer either the forward or the backward 
derivative. One concern is whether our results on the in- 
ferred biases change if we use different derivatives (or a 
2nd order Taylor expansion). To investigate this, we re- 
peated our analysis using backward and "2-sided" deriva- 
tives. In Table |IV[ we show the biases of the three cos- 
mological parameters as, w, £l m inferred from the peak 
counts, with all three types of derivatives. As the table 
shows, the biases from the low peaks remain similar (or at 
least negligible, within errors) for each type of derivative. 
However, high peaks have large and significant differences 
in the biases. 

These difference may indicate a genuine asymmetry in 
the cosmology-dependence of the peak counts, but the 
differences may be partly a numerical artifact due to the 



<78-fi m degeneracy. Given the cosmology-dependence of 
the peak counts (see, e.g., Figure 8 in ref. [24]), it is 
apparent that increases in erg can be mostly compensated 
by decreases in Q m , and vice versa (see also ref. [11]). It 
is easy to imagine that as a result of this degeneracy, 
small changes in the baryonic peak counts that are being 
fit can produce large changes in the as and O m biases - 
depending on where along the as — £l m degeneracy curve 
the baryonic effect forces the best-fit (the residuals could 
then also produce large changes in the w bias). 

In addition to the high peaks, we find that the bi- 
ases from the high-^ 1 power spectrum also change by up 
to 30% between using backward vs. forward derivatives. 
We conclude that, unfortunately, we cannot determine 
the values of the biases of the high peaks and the high-£ 
power spectrum accurately: we need a larger number 
of WL maps to determine the derivatives with better 
precision, and the parameter space needs to be sampled 
more finely with simulations to map out the nonlinear 
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Observable 


derivative 


Aag Aw AQ m 


Low Peaks 


backward 
forward 
2-sides 


0.009 0.01 -0.001 
0.012 0.00 -0.002 
0.009 0.01 


High Peaks 


backward 
forward 
2-sides 


0.061 0.10 -0.024 
0.035 0.05 -0.009 
0.026 0.09 



TABLE IV: The biases of the three cosmological parameters 
as, iv, Q m inferred from the peak counts, using three differ- 
ent types of finite- difference derivatives. The baryonic models 
were again fit by non-baryonic models, using linear interpo- 
lations with either backward, forward, or "two-sided" deriva- 
tives, to make predictions for cosmologies in-between those 
that we simulated. (The biases for fi m in the high peak case 
for the 2-sided derivative could not be determined reliable, ow- 
ing to numerical issues related to the discontinuous switch in 
the derivative type at the fiducial fl m .J 



dependence of the observables on the cosmological pa- 
rameters. Nonetheless, our overall conclusion, namely 
that the low peaks and the law-£ power spectrum are un- 
biased, while high peaks and the high-£ power spectrum 
are noticeably biased, holds for all types of derivative we 
have tried. Furthermore, it would be a remarkable coin- 
cidence if the true biases from the latter two observables 
turned out to be in a degenerate direction; we therefore 
also expect our generic conclusion about the possibility 
of "self-calibration" to remain valid. 

How robust are the results to the number of realiza- 
tions of the baryonic maps? As mentioned above, our 
results are ultimately based on the 200 baryonic maps 
in the fiducial cosmology. To show that these 200 bary- 
onic maps are sufficient to capture the systematic changes 
caused by the baryon cooling, in Table [V] we compare the 
biases of the cosmological parameters inferred from 100 
vs. 200 realizations of baryonic maps. We found that 
the biases from 200 baryonic maps are close to the biases 
from 100 baryonic maps, with the differences of < 0.002 
in Aer 8 , < 0.01 in Aw and < 0.002 in AO m . These dif- 
ferences are smaller than the uncertainties on the biases 
discussed above. This is in reassuring contrast to the ap- 
parently high demand on the accuracy of the derivatives, 
where not even 1000 maps are enough to get stable bi- 
ases for high peaks and the high-^ power spectrum. It 
is also worth noting that when we tried to repeat our 
analysis using only 50 realizations, we found significant 
differences in the biases. This shows that the baryonic 
effects on the peak counts vary from one realization to 
another, and at least >100 realizations are required to 
quantify them reliably. Interestingly, we have found that 
the baryonic effect on the power spectrum is much more 
systematic, with little variation from one map to another 
(see discussion below). 

Effect of noise. An important question is how sensitive 
our results are to the presence of galaxy shape noise. To 
address this issue, we have repeated our analysis with- 
out adding noise to the maps. Naively, one expects that 
noise might suppress the baryon-induced differences in 



Observable 


# of realizations 


Aag Aw AQ m 


Low Peaks 


100 
200 


0.011 0.010 -0.001 
0.009 0.005 -0.001 


High Peaks 


100 
200 


0.059 0.101 -0.023 
0.061 0.103 -0.024 



TABLE V: The biases of the three cosmological parameters 
og, w, Q, m inferred from the peak counts, based on either 100 
or 200 realizations of the baryonic maps. 



the peak counts. This expectation is confirmed by Fig- 
ure [5] which shows, in particular, that the range of peak 
heights that are unaffected is much narrower in the noise- 
less case than in the noisy case. We therefore proceed by 
re-defining the low- and high-peaks in the noiseless maps 
to be those with heights between < K p0 ak < o~o and 
Kpcak > o~Q. Here ctq = 0.022 denotes the standard devi- 
ation of the convergence k in the absence of noise, and 
corresponds roughly to where baryon effects become sig- 
nificant (see the red curve in Fig. [5]). 

The biases of the cosmological parameters inferred 
from the low peaks are found to be (Acs, Aw, AQ m ) — 
(-0.005,0.09,0.001). For high peaks the biases are 
(Acr 8 , Aw, AO m ) = (0.017,0.24,0.006). These noiseless 
bias vectors are generally different from those in the noisy 
case (rows with 200 realizations in Table |v|) . As in the 
noisy case, the biases from the low peaks are still much 
smaller than from the high peaks (we have verified that 
this is also true for the forward and 2-sides derivatives). 
However, the simple expectation that noise reduces the 
bias holds only for both low- and high peaks for w. The 
noise increases the bias in a s for both types of peaks. 
Finally, the f2 m -bias is decreased by noise for low peaks, 
but increased for high peaks. In our analysis, we have 
assumed that the noise is Gaussian, and is known per- 
fectly. The above changes in the bias vectors imply that 
precise measurements of the shape noise will indeed be 
crucial when attempting to model baryon effects in the 
WL data. 

We have also found that the biases in the noiseless 
maps are not as stable with respect to the number of re- 
alizations as in the noisy case. The worst case from back- 
ward derivative is the noiseless low peaks. For example, 
when we use 100 realizations, the biases change by a fac- 
tor of wtwo, to (Acr 8 , Aw, AQ m ) = (-0.011, 0.05, 0.003). 
As we will see in the next paragraph, this is in contrast 
with the behaviour of the bias from the power spectrum, 
which is always very stable to number of realizations. A 
possible explanation of this difference is that changing 
the concentration parameter can impact any individual 
peak either way. This is because the change in the con- 
tribution k to a peak from an individual halo can be pos- 
itive or negative, depending on the impact parameter of 
the halo; halos contributing (especially to the low) peaks 
have a broad range of impact parameters (these points 
will be demonstrated in Figure [9] below) . This effect is 
reduced in the presence of noise, since noise itself tends 
to give the largest contribution of k to low peaks. 
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Robustness of biases from the power spectrum. As men- 
tioned above, we have found that the results from the 
power spectrum, both with and without noise, are al- 
ways very robust with respect to the number of realiza- 
tions of baryonic maps. The small sizes of the error bars 
in Figure [4] already show that baryonic effects have much 
weaker map-to-map variations, compared to those in Fig- 
ure [U The differences in the biases inferred from 200 vs. 
100 realizations are less than 0.001 in Acts, 0.01 in Aw; 
and 0.001 in An, m . In fact, we have found that for the 
power spectrum at £ > 2000, even a single map is suffi- 
cient to determine the biases with <10% accuracy. For 
the power spectrum at I < 2000, which is numerically 
somewhat less stable than the high-£ power spectrum, 
we have found that ^5 maps are needed to achieve the 
same accuracy. These numbers refer to the noisy maps; 
in the noiseless case, the I < 2000 power spectrum would 
require >10 maps. 

The biases from the power spectra are more sensitive 
to derivative types, with absolute changes comparable to 
those for the peaks. However, compared to peaks, the 
biases from the power spectra are larger to begin with, 
making the biases more stable in a fractional sense. The 
differences in biases from the power spectra with and 
without noise are < 30%; the only exception is the noisy 
case and the high-£ range, for which there is a change 
by a factor of 2 in Aw. We have studied the sensitivity 
to the number of bins only for the high-€ case (this is 
because we did not have a sufficient number of pre-bins 
saved in the low and very \ow-£ ranges to increase the 
resolution in those cases) . In our fiducial set of compu- 
tations above, we used 13 high I bins. We found that as 
we increase the number of bins from 13 to 15 to 20, the 
changes in the biases are comparable to those for the high 
peaks, with the 13-to-15-change being a bit smaller than 
the 13-to-20-change. However, again because the power 
spectrum in general has much larger absolute biases, even 
in the worst case, changing from 13 to 20 bins, the cor- 
responding fractional change in the parameter biases is 
low ( <15%). 



B. Cosmology from low-bias statistics 

We have shown that baryon cooling causes a strong 
increase in the number of high peaks and the small scale 
power spectrum. The number of low-amplitude peaks 
and the large-scale power spectrum are both relatively in- 
sensitive to baryon cooling, and thus deliver nearly bias- 
free cosmological constraints. Figure [6] reveals that the 
area of the marginalized 2D error contours from these 
low-biased statistics are roughly comparable to those 
from their more biased counterparts. In this section, we 
will further explore how much of the cosmological sensi- 
tivity resides in the low-bias statistics. 

To assess the fraction of the cosmological sensitivity 
coming from low peaks and large-scale power spectra, we 
first computed a A\ 2 , analogous to the quantity defined 



Observable 


FOM 
(all) 


FOM 
(low-bias) 


Fraction 
(low/all) 


Fraction 
per parameter 


Peaks 


1.30 • 10 b 


0.723 • 10 b 


0.56 


0.82 


Pow. Spec. 


1.10 ■ 10 b 


0.060 • 10 b 


0.055 


0.38 



TABLE VI: The strength of the cosmological constraints for 
peak counts and power spectrum, expressed in terms of a figure 
of merit (FOM; defined as the inverse of the three-dimensional 
error volume in the as, Q. m , W parameter space). FOM values 
are compared for the relatively unbiased low- amplitude peaks 
and very low-l (100 < I < 1000,) power spectrum (second 
column) vs. all peaks and the power spectrum on all scales 
(first column). The third column indicates the ratios of these 
FOM (indicating the fraction of the cosmological constraints 
contained in the low-bias ranges of both statistics), while the 
rightmost column shows the cube root of this ratio (a rough 
estimate of the fraction of the constraint on individual pa- 
rameters). 



in Eq. ^ above, except using the difference AN between 
the mean number of peak in pairs of CDM simulations 
with different Og, w and O m . This measures the signifi- 
cance of the differences caused by the changes in the in- 
dividual cosmological parameters. We find that the A\ 2 
values from the relatively unbiased low-amplitude peaks 
exceed half of the A% 2 values obtained from all peaks. 
On the other hand, the A\ 2 values from the large-scale 
power spectrum are typically reduced by a factor of ~ 5 
compared to using the power spectrum on all scales (the 
exception is w, for which the large angular scales contain 
« half of the total Ax 2 ). 

The above shows that the low peaks contain most of 
the raw cosmological sensitivity to each individual pa- 
rameter, but it neglects degeneracies among parameters, 
which ultimately drive the marginalized constraints. We 
therefore next assess the fraction of the full cosmolog- 
ical sensitivity, coming from low peaks and large-scale 
power spectra, when all three parameters are simultane- 
ously varied. For this purpose, we repeated our anal- 
ysis as outlined in this paper, except we bootstrapped 
the CDM maps from our 45-simulation map set, rather 
than the baryonic maps. This yields constraint in the 
three-dimensional parameter space eg, w and £l m . We 
measure the three-dimensional volumes in the parame- 
ter space containing 68% of the probability for all peaks 
and low peaks separately, and likewise for the very low-£ 
and the full power spectra. The figure of merit (FOM) is 
defined to be the inverse of these volumes, and indicates 
the overall strength of the constraints. The ratio of the 
FOM from low vs. all peaks and the very \ow-£ vs. the 
full power spectra then captures the "fraction" of the full 
cosmological sensitivity in the unbiased statistics. This 
ratio can be converted roughly into a per-parameter es- 
timate by taking the third root, although the differences 
in constraints on individual parameters would of course 
depend on the exact shape of the constraint volume. We 



list these quantities in Table VI 



As is evident from the table, the peaks and the power 
spectrum deliver comparable FOMs when the whole 
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range is used. However, the essentially unbiased low- 
amplitude peaks contain more than half (a fraction 0.56) 
of the FOM, whereas this fraction is only ~ 0.05 for 
the very \ow-£ power spectrum (the latter increasing 
to ~ 0.16 if multipoles up to £ < 2000 are included). 
In terms of the fraction per parameter, the low peaks 
contain approximately 82% of the information, whereas 
the very low-£ power spectrum contains 38% (the lat- 
ter number increasing to 54% with multipoles up to 
I < 2000). We conclude that the low-bias statistics con- 
tain the (slight) majority of the cosmological sensitivity 
for peaks, but a relatively small fraction for the power 
spectrum. 
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C. The insensitivity of low peaks to baryon cooling 

We next examine why the number of low peaks is in- 
sensitive to baryon cooling, as modeled by an increase in 
the halo concentration parameters. In short, we attribute 
the lack of sensitivity to two reasons: halos contributing 
to low peaks have large off-sets from the line-of-sight to- 
ward each peak (and therefore the light rays do not "see" 
the halo cores) and also relatively low masses (so that the 
details of the projected halo density profiles are "washed 
out" by the 1-arcmin smoothing). We demonstrate these 
points explicitly below. 

In Figure [HJ we show the distribution of the impact 
parameters (in units of Riso) a t which light rays, corre- 
sponding to the centers of WL peaks, intersect dark mat- 
ter halos along the line of sight. Following ref. (23], for 
each peak, we identify all halos along the sightline, and 
rank them by their contribution to the peak's height. Go- 
ing down this ranked list (beginning with halos with the 
largest contribution), we sum the k contributions from 
the first n halos until these n halos account for > 50% of 
the total halo contribution to the peak convergence. This 
procedure assigns the number n to each peak (hereafter 
referred to as an "n-halo peak"). 

The top panel relates to high peaks (with heights of 
k > 0.08), which are typically produced by one or two ha- 
los. The panel shows, separately, the impact-parameter 
distribution for the dominant halos for 1-halo peaks, 
and for 2-halo peaks. Low peaks (with heights between 
Cnoisc < Kpoak < 0.08; with a noise = 0.023) are typi- 
cally produced by 4-8 halos. The representative case of 
5-halo peaks is shown in the bottom panel. Our baryonic 
model, with a 50% boost in the concentration parameter, 
was used for this figure (with source galaxies at redshift 
z s = 2 and with galaxy shape noise and 1 arcmin smooth- 
ing included as before). 

As the figure shows, for high peaks, the light rays typ- 
ically go near the halo centers, with impact parameters 
< 0.2i?i8o- In contrast, for low peaks, the light rays 
have large impact parameters, with a maximum near 
~ 0.6i?i8o- This is as expected: as more halos contribute 
to a peak, the contribution from each halo is lower, and 
the halos are less precisely aligned along the line of sight. 
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FIG. 8: The distribution of the impact parameters (in units 
of the virial radius) at which light rays, corresponding to the 
centers of WL peaks, intersect dark matter halos along the 
line of sight. The upper panel is for high peaks, which are 
typically produced by one or two halos. The panel shows, sep- 
arately, the impact-parameter distribution for the dominant 
halos for 1-halo peaks (solid red curve), and for 2-halo peaks 
(dashed green curve), "n-halo peaks" are defined to be those 
for which the sum of n halos along the line of sight accounts 
for at least 50% of the halo contribution to the total peak con- 
vergence. The light rays typically go near the halo centers, 
with impact parameters < 0.2_Rigo- In contrast, for the low 
peaks, shown in the lower panel, the light rays have a large 
impact parameter. These peaks are typically produced by 4-8 
halos; the panel shows the distribution for the typical case of 
5-halo peaks, which has a maximum at ~ 0.6i?i8o- 



Figure [9] shows the fractional difference in the conver- 
gence contribution from halos, caused by a 50% increase 
in their concentration, as a function of the impact param- 
eter (in units of the halo's virial radius Riso)- The three 
curves correspond to halos with masses of 10 12 , 10 13 , and 
lO 14 ft. _1 M0, as labeled. The halos and source galaxies 
are located at redshift z = 0.5 and z s = 2, respectively 
and we assume the halos have NFW profiles, with a con- 
centration parameter (before the 50% boost) given by 
Eq. Q. As the fi gure shows, lower-mass halos are less 
affected by the boost in concentration. These halos are 
more compact, and the 1 arcmin smoothing makes their 
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FIG. 9: The figure shows the fractional difference in the con- 
vergence contribution from halos, caused by a 50% increase 
in their concentration, as a function of the impact parame- 
ter d (in units of the halo's virial radius Riso)- The red plus 
signs, green crosses, and blue snowflakes correspond to halos 
with masses of 10 12 , 10 13 , and 10 14 /i _1 M Q , as labeled. The 
halos and source galaxies are located at redshift z — 0.5 and 
z s = 2, respectively and we assume the halos have NFW pro- 
files, with a concentration parameter (before the 50% boost) 
given by Eq. Q). Lower-mass halos are less affected by the 
boost in concentration. These halos are more compact, and 
the 1 arcmin smoothing makes their convergence contribution 
less sensitive to their intrinsic density profile. 



convergence contribution less sensitive to their intrinsic 
density profile. 

Figures [8] and [9] together allows us to assess the sen- 
sitivity of peak counts to baryons for both high and low 
peaks. High peaks are dominated by halos with higher 
masses 10 13 — 1O 14 /i _1 M , and the dominant halos for the 
l-halo peaks have masses near the upper end of this mass 
range. On the other hand, the low peaks are dominated 
by halos with lower masses, between 10 12 — 10 13 h~ 1 MQ. 
Combining the results from Figures [8] and [9j we see that 
high peaks with masses < 10 14 /i _1 M Q and impact pa- 
rameters < 0.2i?i 80 have typical fractional differences of 
about 8 percent. By comparison, low peaks with masses 
10 12 — 10 13 /i _1 Mq and impact parameters ~ 0.6i?iso, 
have typical fractional differences of about 2 percent, 
much lower than high peaks. This explains why the low 
peaks are less sensitive to baryon cooling. 



V. SUMMARY AND CONCLUSIONS 

In this paper, we have studied the impact of the cool- 
ing and concentration of baryons in dark matter halos 
on WL observables, by manually steepening the den- 
sity profile of each DM halo in a suite of ray-tracing 
N-body simulations. Our main findings are that a) the 
low peaks (0.023 < K pea k < 0.08) are mostly unaffected, 
b) high peaks (K pea k > 0.08) and the power spectrum 
at £ > 1,000 lead to biases if baryonic effects are ne- 



glected in the cosmological parameter estimation, c) the 
bias in high peaks is comparable to the bias for the 
\ow-£ (£ < 2000) power spectrum, and d) the high-^ 
(2000 < £ < 20, 000) power spectrum exhibits a much 
stronger bias in a different direction. 

We find a large increase in the amplitude of the small- 
scale power spectrum, caused by the steepening of halo 
profiles, confirming previous works. It is unsurprising 
that we find a corresponding large increase in the num- 
ber of high peaks, since these peaks are typically caused 
by individual halos. However, the fact that low peaks 
are essentially insensitive to an increase in the halo con- 
centration was surprising (at least to us). This is an es- 
pecially important finding, since these low peaks contain 
the majority of the cosmological information contained 
in the entire set of peaks. We attribute the robustness 
of these low peak to the fact that they are created by an 
alignment of typically 4-8 of low-mass (~ 10 12 — 1O 13 M ) 
halos [2U, with each halo typically offset from the line- 
of-sight towards the peak by a fair fraction of the virial 
radius (~ 0.5i? V j r ). As a result, light-rays corresponding 
to the peak do not pierce the cores of these halos, where 
baryonic effects are most significant. Additionally, the 
halos contributing to the low peaks are compact on the 
sky. As a result, when their intrinsic profiles are con- 
volved with the smoothing filter (with an angular scale 
of 1 arcmin in our case), applied to the raw convergence 
maps, the details of the profiles are washed out. 

We have explicitly computed the biases in the cosmo- 
logical parameters w, as and Vt m when peak counts and 
power spectra are fit neglecting the baryonic effects. We 
find that the biases from the low peaks and from the 
very \ow-£ power spectrum (t < 1000) are small (consis- 
tent with zero within errors) . The biases from high peaks 
and from the low-£ power spectrum (£ < 2000) are signif- 
icant and comparable both in magnitude and sign (e.g. 
both biases are Aw « +0.1 for dark energy). It is an 
important result that the high peaks are no more biased 
than the \ow-£ power spectrum (I < 2000), the latter be- 
ing the range of scales in the baseline plans by the LSST 
WL survey, to infer cosmological parameters [49] . 

Our finding that biases from the high-^ power spec- 
trum (2000 < £ < 2 x 10 4 ) are very large (between 
Aw = —0.3 and —0.6), but in a very different direction 
in the (w, ag, fl m ) parameter space than the bias from 
high peaks, also has an important general implication. It 
suggests the possibility of an effective self-calibration, by 
tuning the concentration parameter to bring the different 
contours in the cosmological parameter space into align- 
ment. However, given the large magnitude of this bias, 
and keeping in mind that additional unknown "baryonic 
parameters" will be present in a fully realistic case, it is 
difficult to see how this self-calibration could be achieved 
to the sub-percent accuracy level targeted by LSST. We 
note, however, that biases from the power spectrum ex- 
tended down to £ < 2, 000 are lower, and still in a mod- 
estly different direction than that of the high peaks, pos- 
sibly allowing more accurate self-calibration. Until bary- 
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onic effects are much better understood, it may be ad- 
vantageous to restrict the analysis to statistics that are 
nearly unbiased, such as the low peaks and the very low 
(£ < 1, 000) power spectrum we identify here. Our result 
suggests than in the case of peaks, less than half of the 
cosmological constraining power may be lost. 

Our two main qualitative results are likely robust: con- 
straints from sufficiently low-amplitude peaks will be rel- 
atively unbiased; and some form of self-calibration will 
likely be realized in a real survey. However, our study is 
clearly highly idealized. We have identified a number of 
caveats arising from our toy model for the impact of as- 
trophysical processes, our limited number of simulations 
(both in terms of sampling the cosmological parameter 
space and in terms of the number of realizations), our 
simplified treatment of source galaxies, noise, and our 
neglect of all instrumental errors. Our results call for 
follow-up work with more extensive simulations and more 
realistic modeling, to address these shortcomings. 
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